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Abstract 

This paper proposes a general framework for probabilistic certification of cancer therapies. The certification is defined in terms 
of two key issnes which are the tumor contraction and the lower admissible bound on the circulating lymphocytes which 
is viewed as indicator of the patient health. The certification is viewed as the ability to guarantee with a predefined high 
probability the success of the therapy over a finite horizon despite of the unavoidable high uncertainties affecting the dynamic 
model that is used to compute the optimal scheduling of drugs injection. The certification paradigm can be viewed as a tool 
for tuning the treatment parameters and protocols as well as for getting a rational use of limited or expensive drugs. The 
proposed framework is illustrated using the specific problem of combined immunotherapy/chemotherapy of cancer. 


1 INTRODUCTION 

The use of dynamic models in the optimization of drug 
scheduling is nowadays a common practice in academic 
works. This long tradition involves different paradigms 
such as optimal control [17,6,12,13,14,2], predictive con¬ 
trol [5], robust control [1] or nonlinear analytic control 
design [10,15]. 

The dynamic models involved in such studies are typi¬ 
cally population models that are built by concatenating 
functional terms (death rate, transition rates, drug ef¬ 
fect terms to cite but few examples). Such models qual¬ 
itatively capture the main phenomena and represent 
their strength and their interaction/coupling through 
dedicated parameters. 

While the qualitative representativity of these models 
is rather easy to assess, the quantitative matching with 
reality strongly depends on the model parameters. The 
latter are unfortunately unknown for a given patient, 
are highly dispersed between patients and vary with 
time and during the therapy for a given patient. 

Some recent works [11,8,1] started attempts to address 
this issue by using robust design in which the therapy is 
computed so that some statement can be obtained for 
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a set of parameters rather than for the single nominal 
parameter vector. A robustness-like statement typically 
takes the following form: 

The scheduled feedback therapy leads to a predefined 
tumor contraction for ANY realization of the vector of 
parameters involved in the model within a predefined 
bounded set 

Therefore, robust design is based on the worst-case 
analysis and can lead to very conservative/pessimistic 
design. This is because the worst case is considered no 
matter how small its probability of occurrence is. 

In order to avoid focusing on few unlikely although very 
bad scenarios, the probabilistic approach seeks state¬ 
ment of the form: 

The scheduled feedback therapy leads to a predefined 
tumor contraction with a probability no less than 
(1 — r])% over all realizations of the parameter vec¬ 
tor assuming that the latter obeys a given probability 
distribution. 

This obviously marginalizes very bad realizations if 
their probability of occurrence is really small. 

This paper formalizes this paradigm for the specific 
case of cancer therapy and gives a complete and under¬ 
standable instance of it in the specific case of combined 
therapy of cancer that involves immunotherapy and 
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chemotherapy. 

It is obvious that given the wide range of problems 
that can be defined in this context, this paper should 
be viewed as an introduction to a rich paradigm and 
a starting point to a large set of variations around the 
necessary specific formulation adopted in the present 
paper. 

The paper is organized as follows: First a general formu¬ 
lation of a class of cancer therapy-related problems is 
given in section 2. Section 3 recalls the framework and 
useful results of randomized optimization approach also 
called the scenario-based approach. The application of 
this framework to the cancer problem defined in Section 
2 is proposed in section 4 in the general case. Finally, 
section 5 fully illustrates the previous sections in the 
particular case of combined immuno/chemotherapy of 
cancer. The paper ends with Section 6 that summarizes 
the paper contribution and gives some hints for future 
investigation. 

2 Probabilistic Certification of a Therapy 

In this section, the concept of a cancer therapy with 
probabilistic certification is clearly stated. 

2.1 The Dynamic Model 

Let us consider a general form of a dynamic system rep¬ 
resenting the evolution of the tumor and the number of 
circulating lymphocytes among other necessary quanti¬ 
ties under a combined action of several drugs injection 
rates u G K”": 

x = F{x,u,p) (1) 

where x G M” is the state of the model while p G 
R"p stands for the vector of parameters involved in the 
model. It is assumed in the remainder of the present pa¬ 
per that 

• Xi stands for the tumor size (to be reduced) 

• X 2 stands for the amount of circulating lymphocytes 
that is commonly used as an indicator of the patient 
health/resistance and therefore, any strategy has to 
be defined such that C(t) > Cmin for all t >0. 

Other state components may be necessary to describe 
the model (namely n > 2) but their exact definition is 
not needed as far as the presentation of the concepts is 
concerned. 

It is assumed that the dynamic model (1) describes the 
evolution of the system under the combined effect of 
different drugs such as chemotherapy, immunotherapy, 
anti-angiogenesis and so on. 



Fig. 1. Temporal structure of the therapy. Example of a 
treatment period consisting of Nt = 3 sub-periods of a duty 
cycle 7 . 

2.2 The Feedback-Based Therapy Protocol 

Let us consider a feedback-based therapy of duration T 
consisting of Nt sub-periods (of duration Ts — T/Nt) 
each of which involving a treatment phase and a rest 
phase as shown in Figure 1 where the injection curves 
have to be interpreted as a multivariable signals when 
several drugs are combined. 


It is assumed that during a treatment period, a sampled 
feedback injection law is used with a sampling period 
T (for instance 2, 4, 6 hours or such) during which the 
injection is maintained constant (see Figure 1): 

u[kT1) = K{x{kT),9c) tS[0,T] (2) 

where x^kr) denote the state of the model at instant kr 
while 9c G K"” is a vector of parameters that are used 
in the definition of the feedback law K. 

In the remainder of the paper, the notation x(k) is used 
instead of x(kT) to simplify the expressions when no 
ambiguity is possible. It is also assumed that the sam¬ 
pling period is a divisor of yTg such that there is an 
integer Ng satisfying: 

iTg = NgT ■ (3) 

It is implicitly assumed that the control law (2) satisfies 
the following saturation constraints: 

K,{x{k),9c)G[Q,ur^[k)] {!,...,n4 (4) 

where Ki stands for the i-th component of K (the z-th 
drug injection value) and where u//^°‘^{k) represents the 
maximum allowable injection rate of the z-th drug during 
the fc-th sampling interval. The fact that the maximum 
injection rate is time-varying is induced by the need to 
meet the constraint on the total amount of available 
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drugs Di for the whole therapy. This constraint can be 
satisfied by using the following definition for 

< min{ui, (5) 

yi{k + l)=y^{k) + TX Ki{x{k),9c) ; yi(0) = 0 (6) 

where yi{k) represents the amount of drug i already 
injected over [0,/ct] meaning that Di — yi{k) is the 
available quantity for the remaining therapy duration 
T — kr. Note that maximum injection rate is also lim¬ 
ited by technical saturation Ui regardless of the amount 
of initially available drug Di. Note that by definition Ui 
is not a design parameter since it is imposed by exoge¬ 
nous technical limitations. 

In section 5, a fully developed example of such control 
law is given for the specific example of combined im¬ 
munotherapy/chemotherapy. For the time being, the 
general non instantiated form ( 2 ) is kept in order to 
preserve the general character of the concepts. 

It comes out that for a given total treatment duration 
T, the therapy is completely defined if the following 
parameters are defined: 

(1) The state feedback parameter vector 9c & Qc, 

(2) The number of sub-periods Nt, 

(3) The duty cycle 7 , 

(4) The allocated drug quantities Di, i e {I,..., n„} 

(5) The tumor contraction ratio 7 c [see (9)] 

These parameters are gathered in the sequel into a single 
decision vector 9, namely: 


Considering a target tumor contraction ratio 7 c at the 
end of the therapy, the success of the therapy can be 
summarized by the fulfillment of the following two con¬ 
straints: 


Xi{T,p,9,xo) 

xi{0) 


< 7c e [o,i[ 


max[Crmn - X2{kT,p,9,xo)] <0 
k 


(9) 

( 10 ) 


as this means that the tumor would be contracted by at 
least 7 c at the end of the therapy while the lymphocytes 
are maintained higher than their lower level Cmin- Ob¬ 
viously, these two constraints can be gathered in a single 
boolean indicator: 


9{0,p) 


0 if (9)-(10) are satisfied 
1 otherwise 


( 11 ) 


where T and xq are supposed to be given and are there¬ 
fore omitted from the list of arguments. Since the indi¬ 
cator g is 1 when the constraints are violated, this indi¬ 
cator is referred to as the failure indicator. 

Definition 2.1 (The Failure Indicator) 

The function g{9,p) defined by (11) is called the failure 
indicator for the therapy defined by 9 and the model de¬ 
fined by the parameter vector p. 

Note however that the constraints (9)-(10) involve the 
UNKNOWN parameter vector p. This makes the definition 
of a successful strategy rather ambiguous. Indeed, two 
statements are possible as mentioned in the introduction 
of this paper: 


9 = (^9c Nt j-fc Di ... Dn^y (7) 

in order to state the probabilistic certification problem 
discussed in the next section. In the sequel, 9 defined by 
(7) is referred to as the therapy design parameter while 
9c is called the feedback design parameter. Therefore, 
the therapy design parameter set includes the control 
parameter 9c but also the time structure and the maxi¬ 
mum injection rates. 


(1) A Robustly certified therapy would be the one 
for which, when using the setting defined by 9, a 
failure indicator g{9,p) = 0 holds for any possible 
realization of p within the admissible set P. 

(2) A (d, ryj-Probabilistically certified therapy 
would be the one for which, when using the setting 
defined by 9 one can state with a probability no 
less than 1 — d « 1 that the expectation of the fail¬ 
ure indicator g{9,p) over P (using the probability 
measure P) is at most equal to ry « 0 . 


2.3 The Concept of Probabilistic Certification 

Note that given the dynamic model (1), the model’s pa¬ 
rameter vector p and the therapy parameter vector 9, the 
evolution of the system can be predicted for any given 
initial state xq (at the beginning of the therapy) so that 
the value of the state x(k) at the sampling instant kr 
can be denoted by: 

x{k) =: X{kT,p,9,xo) ( 8 ) 


For obvious reasons, the parameter d is referred to as 
the confidence parameter (since it is the probability that 
the success statement is wrong) while ry is referred to 
as the precision parameter since it represents the error 
committed w.r.t the ideal achievement g = 0. Note that 
a robustly certified therapy is a ( 0 , 0 )-probabilistically 
certified therapy. 

As mentioned in the introduction, the first concept gen¬ 
erally leads to very pessimistic design since it is based 
on the worst case scenario even if it is very unlikely. 
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Moreover the corresponding computation is extremely 
difficult. The second concept leads to more tractable 
computation and it neglects very unlikely bad scenarios 
leading to more pragmatic design. This is the option 
followed in the remainder of the paper. 


new optimization problem becomes: 


N 

min J (0) under ^9(0, < m 


(15) 


In the above discussion, only the satisfaction of the con¬ 
straints (9)-(10) is considered. As a matter of fact, there 
might be several values of 9 that meet the constraints in 
which case the best 9 should be defined through some 
cost function J{9) to be minimized. In our problem this 
may be 

/ the quantities of the used drugs (which are a part of 
9 according to (7)) 

/ the duty cycle 7 reducing hospitalization periods 
/ any convex combination of the above indicators. 

In the next section, the computational aspect that en¬ 
ables to compute 9 leading to a probabilistic certifica¬ 
tion of the therapy is introduced by recalling the main 
ideas on the general topics of randomized methods [3,4]. 

3 Recalls on Randomized Methods 


which simply replaces the constraint on the probability 
by a different constraint stating that the mean value of 
g{9,p^^^) over N random trials to be lower than m/N, 
or to state it differently that at most m between the 
total number N of trials lead to the violation of the 
specification. It comes therefore that N must be such 
that: 


which is obviously only a necessary condition. This is 
because N must also be sufficiently large so that the 
fulfillment of (15) implies that the condition (14) on 
the probability is satisfied with a probability greater 
than I — (5 with a pre-specified small value 6. that is the 
reason why the minimum value of N that makes this 
implication true involves both the precision specified by 
r] and the confidence specified by 6. 


Consider the following robust optimization problem in 
the decision variable 9 G O C M"'’ and the uncertainty p: 

min J(9) under (Vp) g(9,p)=0 (12) 

where g(9,p) is defined as in (II) by: 

, , ( 0 if specification are satisfied , , 

9(0,P) := / (13) 

( I otherwise 

and where a probability measure P is associated to the 
uncertainty vector p that is assumed to belong to some 
admissible set P. 

The randomized method replaces the original hard 
problem ( 12 ) by the following problem: 

minJ(0) under Pr-p{g(9,p) = 1} < p (14) 


In [3,4], several expressions for the value of N are given 
under different assumptions. In this paper, we are in¬ 
terested in the particular case where the set of design 
parameter 9 is discrete with cardinality uq G N. This 
is because some of the parameters being involved such 
as the available quantities of drugs Hi, the number of 
sub-periods Nt are naturally quantified and cannot be 
viewed as a free real variables. For all the remaining 
variables, one can take some representative values on 
the admissible intervals. By doing so, the optimization 
problem (15) is greatly simplified since it can be solved 
by simple enumeration. Obviously, mixed integer non¬ 
linear programming can also be used following the same 
lines presented in the paper without significant qualita¬ 
tive difference. 

According to [4], in this case, the following proposition 
holds: 

Proposition 3.1 Let m G N be any integer. Let 6 G 
(0,1) be a targeted confidence parameter and p G (0,1) 
be a targeted precision parameter. Take N satisfying: 


where Pr-p{(/(0,p) = 1} represents the probability of 
the event g{9,p) = I (violation of the requirement) 
when p is randomly generated in accordance with the 
probability measure V. 

Now since the computation of the probability term is 
a rather involved and expensive task, the randomized 
method [3,4] simplifies (14) by replacing the probability 
by the mean value over N drawn independent identi¬ 
cally distributed (i.i.d) samples of p in P, namely the 



then any solution to (15) in which the o-t’c ran¬ 

domly i.i.d drawn using the probability measure V satis¬ 
fies the constraint in (If) with a probability > I — d. 

A remarkable property of the expression (17) enabling 
the computation of N is that it is totally independent 
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of the the dimension of p (the number of parameters 
involved in the dynamic model in our application). This 
is of tremendous importance in the context of certified 
therapy since there are typically a quite high number of 
parameters (these are typically the gains associated to 
each functional term in the model). It is a rather good 
news that this does not influence the number of trials 
that is needed to define the constraint in the optimiza¬ 
tion problem (15). This is a rather counter intuitive 
feature for a simple first thought. 

Another interesting feature of Proposition 3.1 is that 
the confidence parameter S appears through a loga- 
rithmique term which means that one can seek highly 
confident assertions without dramatic increase in the 
number of trials. 

In the next section, the use of the randomized method 
summarized in Proposition 3.1 in the certification of 
cancer therapy is presented in the general setting before 
a specific and complete study of a particular case is 
proposed in section 5. 

4 Application to Certification of Therapies 

The application of the framework of the preceding sec¬ 
tion to the certification of cancer therapy can be achieved 
using the following setps: 

(1) Definition of the feedback law: 

First of all, a state feedback law of the form (2) has 
to be designed. This design is problem-dependent 
although some works seek general structures for 
the solution such as in [9]. Another option is to 
adopt systematic use of generic approaches such as 
Model Predictive Control (MPC) [5,16] which can 
be systematically applied as soon as some dynamic 
model (including potentially nonlinear complex 
models), a cost function and a constraint function 
are clearly defined which is the case in our context. 
Note that free open-source available softwares are 
now available for practitioners that enable easy im¬ 
plementation of MPC controllers [7]. nevertheless, 
the definition of the control law delivered by such 
tools still need some parameters to be fixed by the 
user such as the weighting matrices, the constraints 
(total drug available), the sampling period and so 
on. These parameters together with those needed 
to define the time structure of the therapy defined 
in Figure 1 represent what is referred to in the pre¬ 
vious section by the therapy parameter vector 9. 

(2) Definition of the probability measure V 

The feedback law invoked in the preceding item 
accepts the parameter vector p used in the defi¬ 
nition of the model (1) as a given parameter. As 
mentioned in the previous section, the use of the 
randomized method need a probability measure V 


to be defined for use in the generation of the N 
i.i.d set of trials p^^\ £ = 1,..., TV invoked in the 
definition (15) of the relaxed optimization prob¬ 
lem. Three main options are here available: 

(a) In the first, a nominal values can 

be used and the probability measure can be 
defined by a Gaussian distribution around this 
nominal value with a predefined covariance 
matrix. 

(b) Another possibility is to consider that each 
component pi belongs to some interval [p^,Pi] 
and the probability measure represents a sim¬ 
ple uniform distribution (all the values inside 
the interval are treated with equal probability). 

(c) The last option combines the two preceding 
one by adopting Gaussian distributions that 
are saturated by some extreme values p^ and 
pi in order to avoid unrealistic trials to take 
place (such as negative values for a intrinsecly 
positive parameter). 

(3) Definition of the Confidence and precision 
parameters. 

These are the parameters 6 and t] involved in the 
randomized approach. Recall that 1 — <5 defines the 
confidence with which the certification result can 
be assessed while 1 — rj represents the probability 
of failure in the fulfillment of the constraints. Typ¬ 
ical values for these parameters are 6 = 10“^ and 

7 ] = 10 “^. 

(4) Definition of the design parameter set 0. 

This is done by choosing for each component 6i of 
the therapy design parameter a set of representative 
values 

0.:={0f),..., (18) 

covering the presumed interval of interesting values. 
This obviously leads to a discrete set of cardinality: 

ne ng 

ne = Y\_'ni ; 0 = ]^0i (19) 

i=l i=l 

recall that the impact of the value of uq on the 
complexity of the solution (through the number of 
trials TV) appear through a logarithm which means 
that high values of uq can be used to reasonably 
explore the design space. 

We assume that a numbering rule is used inside 0 
so that the uq elements of the discrete set 0 can 
be denoted as follows: 

0 := ( 20 ) 
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ne 

rj = 0.1 

g = 0.05 

q 

o 

II 

g = 0.001 

1 

132 

264 

1317 

13164 

5 

154 

308 

1536 

15354 

10 

163 

326 

1628 

16280 

100 

193 

386 

1930 

19299 

1000 

223 

445 

2225 

22249 

10000 

252 

503 

2515 

25148 


Table 1 

Evolution of the sample size N (number of trials needed to 
achieve the certification) as a function of the precision rj and 
the cardinality ne of the design parameter set 0 (confidence 
parameter S = 10~® is used). 

(5) Computing the sample size N. This can be done 
by choosing an arbitrary value of m (say m = 1) 
and using the above mentioned S, rj and ne in (17) 
to compute N. Table 1 shows the evolution of the 
sample size N (number of trials needed to achieve 
the certification) as a function of the precision ry 
and the cardinality ne of the design parameter 
set 0. The confidence parameter is systematically 
taken equal to <5 = 10“^. 

(6) Draw the model parameter samples. 

Having N at hand, a set of N sample 
i = 1 ,..., is drawn using the probability mea¬ 
sure defined in step (2) above. 

(7) Perform Closed-loop simulations. 

In this step, for each of the ne candidate values 
G 0, cr = 1,..., ne defined in step (4) and 
each of the model parameter vector generated 
in step (6), the resulting model (with used for 
p in (1)) is simulated using the feedback therapy 
defined by This obviously results in • ne 

closed-loop simulation of the model over the ther¬ 
apy duration. Table 1 can be used to evaluate this 
number by multiplying each element of the inside 
matrix (given Af) by the corresponding line value 
of ne. For instance, when using ry = 0.01 and 
ne = 10000, one needs 2515 x 10000 « 25 x 10® 
closed-loop simulations of the therapy. Note how¬ 
ever than with nowadays computers, a single sim¬ 
ulation of commonly used population models takes 
no more than a hundred of microseconds which 
brings the computation time (even for the very de¬ 
manding precision level corresponding to ry = 10“^ 
to less than one hour. 


Remark 1 Note that this estimation of the com¬ 
putational task is pessimistic since the candidate 
values 6^'^) can be visited in a clever way so that 
necessarily unsuccessful values are never tried. For 


instance, if for some quantity of drugs Di and a 
given set of other parameter is unsuccessful, there is 
no need to visit all those combinaison of parameter 
that correspond to lower values of Di. 

Note that for each simulation corresponding to 
p^^'>), the resulting failure indicator: 

:=5(0(-),pW) (21) 

can be computed where g{-,-) is defined by ( 11 ). 
Similarly for any candidate cost function (quantity 
of drugs, duty cycle, etc) the corresponding cost 
matrix can be computed. 

(8) Computing the admissible set of design pa¬ 
rameters 

Having computed the constraints in (15) can 
now be evaluated for each candidate parameter 
by summing the columns of the cr-th line of the ma¬ 
trix namely: 

N N 

= ( 22 ) 

e=i e.=\ 

if the result is lower than m then the candidate value 
0^'^) is considered to be admissible. Therefore the 
admissible set of design parameters is defined by: 

N 

^ := |cr G {l,...,ne} I < mj (23) 

1=1 

(9) Compute the optimal certified therapy 

The optimal therapy is defined by 9^'^ ^ where a* is 
the index of the admissible therapy that minimizes 
the cost function, namely: 

a* = argmin (24) 

(jdA L J 

In the next section, the road map detailed above is ap¬ 
plied to the specific example of combined immunother¬ 
apy/chemotherapy of cancer. 


5 Illnstrative example: Combined immuno/chemo 
therapy of cancer 

The main objective of this example is to enhance a com¬ 
plete understanding of the proposed framework so that 
future works can be initiated using various models, com¬ 
bination of drugs, cost functions and so on. 
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Eq. 

Term 

Description 

(25) 

aa;i(l — bxi) 

Logistic tumor growth 

(25) 

—C 1 X 4 X 1 

Death of tumor due to effector cells 

(25) 

—ksxsxi 

Death of tumor due to chemotherapy 

(26) 

— 5X2 

Death of circulating lymphocytes 

(26) 

— k2X3X2 

Death of lymphocytes due to chemo 

(26) 

S 2 

Constant source of lymphocytes 

(27) 

-70*3 

Exponential decay of chemotherapy 

(28) 

Xi 

9, , X4 

h + xi 

Stimulation of tumor on effector cells 

(28) 

—rx4 

Death of effector cells 

(28) 

-P 0 X 4 X 1 

Inactivation of effector cells by tumor 

(28) 

— k4X4X3 

Death of effector cells due to chemo 


Table 2 

Signification of the terms involved in the dynamic model 
(25)-(28) [source [10]] 

5.1 The dynamic model 

Consider the dynamic model used in [10] in which a ex¬ 
ternal source of effect or-immune cells can be adminis¬ 
tered in addition to the chemotherapy drugs. The model 
involves 7 states and 2 control inputs that are defined as 
follows: 


5.2 Definition of the feedback law 


The starting point in the design of the feedback law lies 
in the fact that according to (27) if one can guarantee 
that X 3 always satisfies the inequality 


X 3 > X 


max 

3 


:= max 


{x,l3) 

/q l^2ildCmin 


— X 2 ) + 5 x 2 — S 2 
-k 2 X 2 


(32) 


for some /3 > 1 then according to (26), the evolution of 
the lymphocytes population would satisfy the inequality: 

X2 > M2(/3C'™n - X 2 ) (33) 

which simply would imply that as soon as X 2 becomes 
lower than pCmin > Cmin then it can only increase. 
This obviously prevent the health constraint X 2 > Cmin 
from being violated. 

The next step is to observe that meeting the inequality 
(32) on X 3 can be guaranteed if one uses equation (27) to 
induce a corresponding limitation in the chemotherapy 
drug delivery. This can be done by if the following con¬ 
straint is satisfied on the chemotherapy drug injection 
rate U 2 '. 

U2<lox'f'°‘^{x,P) (34) 


xi tumor cell population 

X 2 circulating lymphocytes population 

X 3 chemotherapy drug concentration 
Xj^ effector immune cell population 
X 3 quantity of already delivered chemo drug 
xq quantity of already delivered immuno drug 
xy remaining time for therapy 
Ml rate of introduction of immune cells 
U 2 rate of introduction of chemotherapy 

The dynamic model takes the standard form (1): 


xi = axi(\ — bxi) — C1X4X1 — k^x^xi ( 25 ) 

±2 = —5X2 — k2X3X2 + S2 ( 26 ) 

X3 = -I0X3 + U2 ( 27 ) 

Xi 

±4 = g- - X4 — rx4 — P0X4X1 — kiX4X3 + SiUi ( 28 ) 

h + xi 

X 5 =ui ; 0:5(0) =0 ( 29 ) 

xq=U 2 ; 0:5(0) = 0 ( 30 ) 

±7 =-I ; 0:5(0) =T ( 31 ) 


This constraint has to be combined with the other con¬ 
straints ( 6 ) imposed on ui in order to meet the technical 
constraint ui < ui and the one induced by the limited 
amount of chemotherapy drug that is available for the 
therapy. This leads to the following definition of : 

ur^ix,P) = min{u2,joxr^ix,f3), ( 35 ) 

I 700:7 J 

where the last term comes from ( 6 ) in which, the already 
injected chemo drug 0:5 and the remaining time for the 
therapy 0:7 are used. As for the immunotherapy drug, 
the following simple definition is used since no other 
limitations are to be considered: 

(36) 

700:7 J 

From now on, when we write u = ^““^( 0 :,/3), this is 
to be interpreted component-wise using (35)-(36). Note 
that the above definitions involve the first parameter 
/3 > 1 that is a part of the control design parameter 9c. 


Ml, 


where the description of the role of each groups of term 
is given in Table 2. Note that the dynamic model (25)- 
(28) involves Up = 14 parameters. The nominal values 
of these parameters as used in [ 10 ] are summarized in 
Table 3. 


The bounds defined above gives the maximum values 
that are possible to be administered. The effectively 
applied values are defined according to the targeted 
tumor decrease. More precisely, assume that an expo¬ 
nential decrease is targeted. Such decrease would be 
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param 

value 

param 

value 

param 

value 

a 

4.31 x 10"® day-^ 

b 

1.02 x 10“^"* ceir^ 

Cl 

3.41 x 10“^® [eell ■ day)~^ 

/ 

4.12 x 10"^ day-^ 

9 

1.5 x 10"^ day-^ 

h 

2.02 x 10^ cell‘d 

k2, ka 

6 x 10~^ day~^ 

fci 

1 

1 

0 

X 

00 

Po 

2 x 10“^^ [eell ■ day)~^ 

Si 

1.2 x 10"^ cell ■ day~^ 

S2 

7.5 x 10® cell ■ day ^ 

S 

1.2 x 10"^ day-^ 

7 

1 

7 

0 

X 

0 






Table 3 

Nominal values of the 14 parameters involved in the dynamic model (25)-(31). 


characterized by the following condition: 


stopped, otherwise, 



Xi 


(37) 


Now for a continuous decrease, when r = 3/T one can 
achieve a settling time (at 0.05 of the initial value) at 
the end of the therapy. In the proposed therapy proto¬ 
col however, because of the potential rest period, much 
higher value of r would be necessary to achieve the same 
contraction and this value strongly depends on the drug 
delivery periods Tg and the duty cycle 7 . That is the 
reason why r is supposed to be the second component 
in the control design vector Oc- 

Denoting by -Fi(x) the r.h.s of (25), the ideal condition 
(37) becomes: 

E(x,r) := < —r (38) 

Xi 


The idea is then to define the feedback using the hys¬ 
teresis that is defined in terms of the function E{x,r) 
as shown in Figure 2. More precisely, the feedback is de¬ 
fined in terms of E{x(k), r) and its past value (over the 
past sampling period, namely E{x{k — l),r) by: 


If xi > threshold 

' ^max II E> -ar 
0 if f; < -r 

u := < 

^max if ^ —ar) and AE < 0 

0 if FI e (—r, —ar) and AE > 0 
else u = 0 


(39) 

(40) 

(41) 


where AEik) := E{x{k),r) — E{x{k — l),r) and where 
a € ( 0 , 1 ) is a design parameter (this is the third pa¬ 
rameter defined so far as a component of the control 
design parameter vector 


The rational behind this definition can be understood 
based on the following comments: 


/ If FI > —ar, this is interpreted as the tumor is not 
decreasing enough, then the maximum drug intensity 
is used. 

/ If FI < —r, this is interpreted as the tumor is decreas¬ 
ing fast enough and the drug delivery is interrupted 
to privilege the patient health and to save drugs. 

/ The remaining hysteresis-like rule are used to define 
the control level over (—r, —ar) 

To summarize the discussion regarding the definition of 
the feedback law, it comes out that the vector of control 
design parameter can be defined by: 

9c := (/3 r a)^ G [l,oo] x [0, 00 ] x (0,1) C (42) 

and more realistic set 0 c to which these parameters can 
be defined by: 

0c := {1.05,2} X {0.05,0.25,0.5,0.8} x {0.1,0.5,0.8} 

(43) 


which is a set of cardinality 24. 

To this set of options, we have to add the extra pa¬ 
rameters Nt, 7 , 7 c, Di and D 2 that are involved in 
the definition of the therapy parameter 9 [see (7)]. For 
this example, the current specific choices will be used 
regarding these design parameters: 

/ The contraction factor 7 c = 0.1 is fixed. 

/ Duty cycle 7 G {0.3,0.5,0.8} 

/ Number of sub-periods = {4, 6 } 

/ Available quantities of drugs: This is parametrized to 
be a quantized fraction of the maximum injectable 
quantity given the therapy duration T and the duty 
cycle 7 : 

Di:=dxjTu^ dG {0.25,0.5,0.75,1} (44) 


/ If the tumor is too small then the treatment is 


where 'jTui is the total amount of drug that is possible 





u = K{x] 



yTTlaX 











->- 


- > 


-r -ar 0 E{x,r) 

Fig. 2. The definition of the feedback law. Note that 
^max ^ 1^2 that the curves have to be interpreted compo¬ 
nent-wise where is given by (35)-(36). 

to inject given T, 7 and the maximal intensity bounds 

Ui. 


This new set of parameters is of cardinality 24 which 
together with the set of control parameter leads to a 
total cardinality 


ne = 24 X 24 = 576 

In the sequel, the minimum number of circulating lym¬ 
phocytes is taken equal to Cmin = 5 x 10^. The bounds 
Ml = 50 and M 2 = 1 involved in (35) and (36) are 
used. The total duration of the therapy is taken equal 
to r = 60 days. The threshold involved in (39) below 
of which the treatment the drug injection is stopped is 
fixed in the sequel to 10“^. The sampling period used to 
update the feedback is taken equal to r = 4 hours. 

Before getting into the certification issue, Figure 3 
shows the results of the therapy using different set of 
therapy design parameters. These plots show how the 
choice of the design parameters systematically meet 
the constraint on the minimum level of circulating lym¬ 
phocytes as it should be expected from the control law 
design while the contraction of the tumor and its inten¬ 
sity strongly depend on the parameter choices. All the 
scenarios start from the common initial state: 

xo = (5 X 10®, 10®, 0,10®, 0, 0, t) 


5.3 Generation of the model parameters sample 

Using the confidence level defined by <5 = 10“®, the 
precision level defined by d = 0.01 and the cardinality 
Me = 576 to compute the sample size N using the ex¬ 
pression (17), it comes out that the sample size if given 
by 

N = 2155 



Simulation 1. Results with the design parameters r = 0.5, 


A = 0 . 757 rMi, 7 = 0.4, /3 = 



1.05, Nt = 6 and a = 0.5 



Simulation 2. Results with the same parameters as in Sim¬ 
ulation 1 but with the duty cycle 7 = 0.7 instead of 0.4. 

Tumor (cell) Circ. lymphocytes (cell) 

• 10 ® - 10 ® 



Simulation 3. Results with the same parameters as in Sim¬ 
ulation 2 but with the hysteresis parameter a = 0.1 instead 
of 0.5. 



Simulation 4. Results with the same parameters as in Simu¬ 
lation 2 but with the available drug quantities Di = 0.5yTui 
instead of A = O.Tb'yTui 


Fig. 3. Four different simulated therapies with four different 
sets of therapy design parameter vector 0. 


Given the value of uq = 576, this means that the com- 
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putation of the best therapy design parameter 6 needs 

Number of simulations = 1,241, 280 

and since a single simulation of the feedback therapy 
over the therapy duration T = 60 days costs approxima- 
tively 70 /r sec, it comes out that the whole computation 
of the parameters for a certihed therapy can be done in 
approximatively 90 seconds. 

Regarding the definition of the probability measure, as 
mentioned in section 4, there are several ways to define 
how the true parameters spread around the nominal val¬ 
ues given in Table 3. The one used hereafter to illustrate 
the methodology considers that each parameter pi of 
the model has a uniform probability over the following 
interval that includes the nominal value 

Pi € [Ai, A 2 ] X (Ai, A 2 ) G (0,1) X (1, 00 ) (45) 

More precisely, if the pair Ai = 0.6 and A 2 = 1.8 are 
used, this means that the probabilistic certification holds 
when each parameter can take with equal probability 
any value in the interval between 60% the nominal value 
and 180% of the nominal value. 

5.4 Validation 

Several validation scenarios are proposed in this section 
depending on: 

(1) The level of uncertainties: Three couples of (Ai, A 2 ) 
are used leading to three uncertainty levels: 
[-10%,-bl0%], [-20%,-b20%] and [-40%,-b80%] 
which correspond to the pair (Ai,A 2 ) given by: 
(0.9,1.1), (0.8,1.2) and (0.6,1.8) respectively. 

(2) The criterion that is used to define the optimal pa¬ 
rameter over the admissible set of values, namely, 
two criteria are used: 

(a) the minimization of the quantity of drugs. This is 
done by minimizing the parameter d involved in 
the definition (44) of the quantity of drug avail¬ 
able for the whole therapy. 

(b) The minimization of the hospitalization periods. 
This is done by minimizing the parameter 7 which 
is the fraction of the treatment according to Fig¬ 
ure 1 . 

The objective is to show how the optimal therapy pa¬ 
rameters are affected by these above paradigms leading 
to different but certified therapies over all possible real¬ 
izations of the model’s parameters. 


Uncertainties 

Min drug 

Min Hospitalization 


q 

II 


q 

II 



r = 0.25 


r = 0.5 


[- 10 %,-bl 0 %] 

a = 0.5 

7 = 0.8 


a = 0.8 

7 = 0.3 



II 


II 



= 0.25; 


\d = Q.7^) 



q 

II 


II 



r = 0.25 


r = 0.5 


[- 20 %, -b 20 %] 

a = 0.5 

7 = 0.8 


a = 0.5 

7 = 0.3 



II 


II 



\ d = 0.5 / 


Vd = 0.75/ 



II 


II 



r = 0.05 


r = 0.25 


[-40%, -b80%] 

a = 0.5 

7 = 0.8 


a = 0.5 

7 = 0.5 



II 


II 



\d = 0 . 5 / 


\d^i ) 



Table 4 

Optimal therapy design for different level of model uncer¬ 
tainties and different cost function. 

(1) Higher uncertainties implies higher values of /? as 
this parameter is used in (32) to consider that the 
lower bound is pCmin rather than the real lower 
bound Cmin- In that sense, /? allows for a security 
margin on the constraint satisfaction. 

(2) Minimizing drug and minimization hospitaliza¬ 
tion seem to have opposite effects, at least for the 
feedback design considered in the present paper. 
Indeed, higher values of 7 enable to reduce the os¬ 
cillation in the tumor size that are induced by high 
periods of drug-free rest and hence use less total 
amount of drugs. 

(3) This last comment also explain why in the presence 
of high uncertainties, higher number of sub-periods 
becomes mandatory in order to reduce the drug- 
free rest periods. 

Note that some of these comments probably hold only 
for the feedback strategy adopted in the paper which 
is simply given here for the sake of illustration of the 
general certification methodology. 


Table 4 shows the optimal therapy design for these six 
different contexts. Several comments may help for a bet¬ 
ter understanding of the results shown in this table: 
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Finally, Figure 4 shows the validation of the certified 
strategies over a sample of scenario containing 5 times 



Uncertainty [-40%,+80%], Min Drng 



Uncertainty [-40%,+80%], , Min Hospitalization 



■> -Di/x5(r) I 

. D2lx^(T) j , 

0 mint^T[X2{t)/C,nm] \ 

— 7c i 

— 1 _ j 

10 '“ 10 * 10 * 10 * 10 * 10 ° 10 * 10 * 
Tumor contraction xi (T) jx\ (0) 


(a) Uncertainty [—40%, +80%] 


Uncertainty [-20%,+20%], Min Drug 



o Dyjx^^T) 
• D2/xi,{T) 


minteT[x2{t)jC,nin] \ 

-7c I 

- 1 I 

10 * 10 ® 10* 10 * 10 * 10 ' io“ 

Tumor contraction a;i(r)/a;i(0) 


Uncertainty [-20%,+20%], , Min Hospitalization 



c Dilxi{T) 

. D2lX6{T) 

■■ minteT[x2{t)/Cm 
-—7c 


Tumor contraction a;i (T) /x\ (0) 


(b) Uncertainty [—20%, +20%] 


Fig. 4. Validation of the certified optimal therapies over 5 x V « 11, 000 scenarios. The fact that almost all the dots appears 
in the upper left corner means that the contraction level = 0.1) is achieved, the health constraint is satisfied and that no 
more than the allowable is used. 
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more scenarios than those used in the sample of size 
N for the optimization purposes. This corresponds to 

5 X 2155 = 10775 scenarios. The fact that almost all 
the dots belongs to the upper-left corner means that 
the tumor contraction by at least a factor of 7 c = 0.1 is 
achieved, the health constraint is satisfied and that the 
quantities of drug used during the therapy is lower than 
the allowable one. the fact that sometimes the quantity 
of drug used is 10 times smaller than the available one 
comes from a specific combination of parameters in the 
interval that makes the decrease of the tumor possible 
without much drug injection. 

6 Conclusions and Future Works 

In this paper, a general framework is proposed for the 
probabilistic certification of combined therapy of can¬ 
cer under tumor contraction and health constraint. The 
proposed solution is based on the randomized method 
that enables to transform the standard robust worst- 
case approach by a tractable problem with probabilistic 
constraints. The general concepts introduced are illus¬ 
trated in the specific case of combined immunother¬ 
apy/chemotherapy of cancer. 

The framework proposed in this paper can be used ei¬ 
ther to define the level of confidence that can be affected 
to a therapy with a given protocol and a given available 
quantity of drugs; or to determine what is the quan¬ 
tities of drugs and what is the protocol to be used in 
order to achieve a targeted level of confidence. As such, 
the framework can be viewed as a decision making tool 
that enables different options to be compared based on 
reliable computation. 

As mentioned earlier, this contribution can be viewed 
as a starting point for a completely new approach to 
model-based control of tumors since it compensates for 
the oversimplifying character of population models by 
allowing high level of uncertainty on the value of the 
model’s parameters. 
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